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Steady solitary and generalized solitary waves of a 
two-fluid problem where the upper layer is under 
a flexible elastic sheet are considered as a model 
for internal waves under an ice-covered ocean. The 
fluid consists of two layers of constant densities, 
separated by an interface. The elastic sheet resists 
bending forces and is mathematically described by 
a fully nonlinear thin shell model. Fully localized 
solitary waves are computed via a boundary integral 
method. Progression along the various branches of 
solutions shows that barotropic (i.e. surface modes) 
wave-packet solitary wave branches end with the 
free surface approaching the interface. On the other 
hand, the limiting configurations of long baroclinic 
(i.e. internal) solitary waves are characterized by 
an infinite broadening in the horizontal direction. 
Baroclinic wave-packet modes also exist for a large 
range of amplitudes and generalized solitary waves 
are computed in a case of a long internal mode 
in resonance with surface modes. In contrast to the 
pure gravity case (i.e without an elastic cover), these 
generalized solitary waves exhibit new Wilton-ripple- 
like periodic trains in the far field. 



1. Introduction 

Internal gravity waves are ubiquitous in the ocean 
owing to density stratification arising from salinity and 
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temperature variations. Waves are often generated by tidal or other currents flowing over 
topography in the stratified ocean. Internal waves play an important role in transferring heat, 
energy and momentum in both horizontal and vertical directions and, owing to relatively small 
density differences in the fluid, are much lower than surface waves at the interface between 
water and air. Internal interfacial waves arise in situations when the stratification has sharp 
density variations such as at the lower boundary of the surface mixed layer of the ocean, 
at the interface between fresh and sea water at the mouths of rivers, or in exchange flow 
situations such as in the Strait of Gibraltar. A simple mathematical idealization for studying 
internal interfacial waves is the wave propagation on the sharp density discontinuity between 
two immiscible fluids. The 'rigid lid' approximation is commonly used for this two-layer system, 
which means the top of the upper layer is bounded by a horizontal rigid wall. However, the 
well-known 'dead-water' phenomenon, first reported by Fridtjof Nansen whose ship experienced 
an unexpected strong resistance when he was sailing on calm seas in 1893, implies that the 
interaction between the upper free surface and the internal interface in a density-stratified fluid 
may be important. Ekman [1] produced a series of elegant laboratory experiments on the dead- 
water problem, and explained that the force exerted by the ship to propel itself forward was 
instead translated into the generation of long internal waves which caused an apparent drag 
on the ship. 

In the polar region, floating sea ice is subjected to wave motion and can both generate and 
respond to the internal waves. The wind-driven small ice pack is a source of internal wave 
generation in the form of pressure ridge keels (see [2] and the references therein). A heavy 
mass moving on a large floating ice sheet with small velocity can excite internal waves if the 
underlying water is stratified, and may experience an anomalous drag analogous to the dead- 
water phenomenon (see [3] for details). Furthermore, the ice cover, which can suppress the 
small-scale surface wave noise such as wave turbulence and wind ripples, can provide a surface 
signature for internal waves. Czipott et ah [4] reported coherent measurements of ice tilt forced 
by internal waves in the Arctic Ocean. This paper is motivated by the physical significance of 
the coupling between internal waves and sea-ice cover in the polar region. This work is devoted 
to the numerical study of the permanent coherent structures — solitary waves and generalized 
solitary waves — in an idealization for the physical problem, namely a two-layer fluid system with 
a density interface and covered by a flexible elastic sheet that resists flexural motion. 

Interfacial waves have been intensively investigated over the past several decades in 
modelling, theory and numerics. Under the 'rigid lid' approximation and the long wave 
assumption (the thickness of one or both layers are small compared with the characteristic 
wavelength), different model equations can be derived to study two-layer interfacial gravity 
waves. On the numerical side, for travelling waves, Meiron & Saffman [5] and Turner & Vanden- 
Broeck [6] first computed fully nonlinear two-dimensional gravity interfacial waves, including 
both periodic and solitary waves. 

More recently, some three-dimensional studies have included either interfacial effects or 
flexural surface effects. For example, fully localized three-dimensional interfacial gravity- 
capillary solitary waves were found by Parau et ah [7]. Three-dimensional surface flexural gravity 
waves are considered, in the presence of a pressure forcing in Parau & Vanden-Broeck [8], in the 
asymptotic limit of wave-packet dynamics in Milewski & Wang [9] and localized solitary waves 
are described in Wang et al. [10]. 

There are relatively fewer studies on the interfacial waves with free surface. For the pure 
gravity waves, Moni & King [11] computed barotropic solitary waves — waves in which both the 
interface and free surface are elevation waves — using the irrotational Euler equation. Michallet & 
Dias [12] and Parau & Dias [13] computed solitary waves with non-decaying periodic tails which 
are usually called generalized solitary waves. An extensive numerical study on interfacial solitary 
waves with free surface was recently carried out by Woolfenden & Parau [14] when the capillary 
effect was included either on the free surface of the upper layer or on both the top surface and the 
interface. Theoretical work by Iooss [15] has considered the spectrum of the linearized operator 
near equilibrium when the lower layer is infinitely deep. A number of model equations have been 



introduced and investigated in the spirit of long wave or Boussinesq approximation, including for 
example Dias & Il'chev [16], Fochesato et al. [17] and Craig et al. [18]. 

In this paper, we numerically study steady interfacial gravity waves in a two-layer fluid system 
with an upper free surface consisting of a flexible elastic sheet modelling the floating ice cover. 
For the sheet, we use a fully nonlinear elasticity model based on the special Cosserat theory of 
hyperelastic shells proposed by Toland [19] for the study of flexural-gravity waves. Therefore, 
there are two restoring forces across the elastic sheet: gravity and an elastic bending force which 
appears as a pressure jump in the Bernoulli equation at the top surface. 

The paper is structured as follows. In §2, we briefly present the mathematical formulation of 
the full problem and its linear behaviour. In §3, we describe the boundary integral method and the 
series truncation method, and thereafter present numerical results for steady symmetric waves, 
including both solitary waves and generalized solitary waves. In the conclusion, we summarize 
our new findings and discuss possible extensions to the work. 



2. Formulation 

(a) Mathematical formulation 

Consider a two-dimensional incompressible, inviscid fluid bounded below by a rigid horizontal 
wall and above by an elastic sheet. For the sake of simplicity, we assume that the fluid is composed 
of two immiscible layers separated by a sharp interface. The density in each layer is assumed to 
be constant: a lighter fluid with density p + lies on a heavier fluid with density p~ (figure 1). 
Cartesian coordinates (x,y) are introduced such that the y-axis points vertically upwards, and 
the x-axis lies in the line of the undisturbed interface. In each layer, the flow is irrotational, and 
the deformations of the elastic cover and of the interface between the two layers propagate at 
a constant velocity c. Therefore, we can choose a frame of reference moving with velocity c and 
reduce the problem to a steady one. Let h + and h~ be the undisturbed depths of the upper and 
lower layers, respectively, and denote the equations of the perturbations of the interface and of 
the upper free surface by y = rj~(x) and y = r) + {x) + h + . 

A natural non-dimensionalization of the whole system involves the depth of the upper 
layer h + and the speed c as reference quantities. Therefore, we define the following non- 
dimensional numbers: 

H = h — l R= — <1 and F= — =, (2.1) 



where g is the acceleration of gravity and F is the Froude number. For the sake of simplicity, we 
continue denoting by 77 ± the disturbance of the top free surface and the interface after rescaling. 
The velocity potentials 0 + in the upper layer and 0~ in the lower layer are both governed by the 
Laplace equations 

+ 0yy = 0^ -H < y < yf (x) (2.2) 

and 

<t>xx + 4>yy = 0, rj-(x)<y<l + rj+(x). (2.3) 

At the rigid bottom y = —H, the slip boundary condition reads 

4>-=0. (2.4) 

There are two kinematic boundary conditions and a pressure balance condition on the interface 
y = rj~(x), which are given by 

= ( t>x 1 lx' ^y = ^tVx ( 2 - 5 ) 

and 

F 2 

— [|V0- 1 2 - R\ V0+| 2 - (1 - R)] + (1 - R)n~ = B- (2.6) 
where B~ is a Bernoulli constant. 




y = -H- 



Figure 1. Schematic description of the two-layer fluid geometry after non-dimensionalization. 



Finally for the top free surface, the two restoring forces are gravity and the flexural elasticity 
of the elastic cover. We neglect the inertia and the thickness of the plate. These two effects can be 
easily incorporated into the system for steady problems (see [20] for the expression of the inertia 
and [21] for how to model the system with the thickness of the plate). The pressure jump exerted 
by the ice sheet owing to flexing — even with these approximations — has been modelled with a 
variety of methods (see [9] for a comparison). In this paper, we use the Cosserat model introduced 
by Toland for this problem (see [19] for the detailed derivation). In contrast to the Kirchhoff- 
Love model, this model supports both elevation and depression solitary waves. Therefore, the 
kinematic and dynamic boundary conditions at the free surface y = 1 + ri + {x) take the form 



Px 



(2.7) 



and 



y [|V0+| 2 - 1] + n + + Eb ( dss* + y 



= B+. 



(2-8) 



Here, B + is another Bernoulli constant, k is the curvature of the free surface deformation, s is the 
arc-length parameter and 



Eh 3 



12(1 - v 2 )p+g(h-+ 



is the dimensionless flexural rigidity of the ice sheet, where E is the Young modulus, v is the 
Poisson ratio and h is the thickness of the elastic sheet. 



(b) Linear theory 

We linearize the whole system (2.2)-(2.8) around the trivial uniform stream solution: B~ = B + = 0, 
(p~ = 0 + = x and rj~ = r) + = 0, i.e. for arbitrary positive wavenumber k, we substitute the regular 
expansion ansatz 

(p~ = x + €(j)~ cosh[k(H + y)] sin(/cx) 

4> + = x + 6[0^ cosh(/cy) + sinh(fcy)] sin(/cx) 

rf~ =efj~ cos(kx) 

r) + = efj + cos(kx) 

into the equations (2.5)-(2.8) and drop the terms of 0(e 2 ) and higher. We then obtain a 
homogeneous linear algebraic system for 0 _ , fj~ and By requiring the solution for 



this system to be non-trivial, we obtain, after some algebra, the linear dispersion relation in the 
form of 

b(k) ± Jb 2 {k) - ±a{k)c{k) 



with 



and 



(F ± ) 2 - ■ 

v ) 2ka(k) 
a(k) = 1 + R tanh k tanh(fcH), 

b(k) = tanh/c + tanh(fcff) + E b /c 4 [tanh/c + Rtanh(kH)] 
c(k) = (1 - R)(l + E h k 4 ) tanh/c tanh(fcff). 



(2.9) 

(2.10) 
(2.11) 
(2.12) 



There are two distinct Froude numbers F + and F~ in the dispersion relation which are referred 
to as the external mode (also called fast mode) and internal mode (also called slow mode), 
respectively. It is obvious that (F + ) 2 -> oo and (F~) 2 -> 0 as k -> oo. 

We first consider the external and internal modes in the long wave regime. Then, the 
wavenumber k is close to 0, and these two modes can be expanded as 



and 



where 



(F+) 2 = 1 (1 + H + Ax) - a + k 2 + 0(k 4 ) 
(F~) 2 = j(l+H — Ax) — a~k 2 + 0(k% 



RH(1 + H + ^ 1 ) + 1 ^ + ^ 



RH(1+H-Ai) + 



1+H 3 



Ax 



A x = ^(l + H) 2 - 4H(1 - R) 



and 



A 2 - 



(H - 1)(H 3 - 1) 



+ 2RH 



1+H 2 



+ (1 - R)H 



(2.13) 
(2.14) 

(2.15) 

(2.16) 

(2.17) 
(2.18) 



It follows that a + is positive definite for arbitrary 0 < R < 1, but the sign of a - is negotiable and 
depends on the choices of H and R (it is not difficult to verify this by considering the case of 
large H). This fact implies that the minimum of the external mode always occurs at a non-zero 
wavenumber (typical dispersion curves for F + are presented in figure 2a). In figure 2b, we present 
three dispersion curves for internal modes. It shows that the maximum of the internal mode 
can occur at a non-zero wavenumber or at zero wavenumber depending on the choices of the 
parameters H, R and Eb- 

Another terminology also commonly used in the theory of multi-layer fluids is that of 
'barotropic' and / baroclinic / modes. In the case of a two-layer fluid with a free surface, the 
barotropic mode is the one whose waves on the free surface and on the interface move in phase 
with each other (i.e. both interfaces are elevated or depressed simultaneously). On the other hand, 
if the two waves have opposite phases, they are denoted as the baroclinic mode. In linear theory, 
the phase relationship between the top free surface and the interface is determined by the sign of 
y) + /fj~. This can be expressed as 



: cosh(/c) — 



[1 - R - kcoth(kH)F 2 ] sinh(fc) 
kRF 2 ' 



We now consider this ratio in the long wave limit. This gives 



lim — 

k^0+ 7]~ 



F=F± 



RH(F 



±\2 



2(1 + RH)H(1 - R) 



(1+H)=f V(1+H) 2 -4H(1-R) 



■H(l-R) 



(2.19) 



(2.20) 
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Figure 2. (a) External (fast) mode for R = 0.9, £ b = 0.5. Solid line: H = 2; dashed line: H = V, dotted line: H = 0.5. 
(b) Internal (slow) mode for R = 0.1, H = 2. Solid line: £ b = 2.5; dashed line: £ b = 1.4755; dotted line: £ b = 0.5. 
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Figure 3. The nature of the modes at the extrema (described in the text) of F ± (k). (a) Fast modes F + with £ b = 0.5. Solid line: 
H = 3; dashed line: H = | . The large graph shows that the modes are barotropic, whereas the inset shows the extrema occur 
at finite k. (b) Slow modes F~ . Solid line: H = \, £ b = 0.1; circles: H = 3, £ b = 2. The large graph shows that the modes are 
baroclinic, whereas the inset shows the extrema switch between k = 0 and finite k. 



By using (2.13), (2.14), (2.17) and (2.18), it is easy to show that the right-hand side of (2.20) is 
positive for F + and is negative for F~ . It follows that the internal modes are baroclinic and the 
external modes are barotropic in the long wave regime. 

We now consider the regime for wavenumbers away from k = 0. In terms of localized travelling 
waves, the most interesting points on the dispersion curves are the global minimum of the 
external mode and the global maximum of the internal mode. Gap solitary waves are likely to 
be found between these two extrema, otherwise they would resonate with linear periodic waves, 
resulting generically in non-decaying tails in the far field. The value of k at which these extrema 
are attained is called 'critical k' . Figure 3 illustrates the behaviour of the critical k and the types of 
waves that occur (baroclinic or barotropic) for some values of the physical parameters. By fixing 
E b = 0.5 and varying R between 0.2 and 1, figure 3a shows that the minimum of F + is always 
attained at non-zero wavenumber (see the inset in figure 3a where we plot the critical value of 
k versus R), and that the external modes are barotropic for both H = 3 (solid line) and H = ^ 



(dashed line). Figure 3b considers the maximum of F~ . For a small coefficient of flexural rigidity, 
it occurs at zero wavenumber, and the quantity of (2.19) is negative-definite from (2.20) (see the 
solid line in figure 3b with E\> = 0.1 and H = ^). For another set of parameters (E\> = 2 and H = 3), 
there is a jump at R ~ 0.15 where the maximum of F~ transits from a non-zero wavenumber to 
zero wavenumber as the density ratio increases, whereas the modes at these maxima are still of 
baroclinic type (see the circles in figure 3b). 



3. Numerical results 



(a) Solitary waves 

(i) Boundary integral method 

Here, we describe the numerical method for the fully localized travelling wave solution of the 
nonlinear problem defined by (2.2)-(2.8). The hodograph transform is used in the numerical 
scheme, and the auxiliary function (f below) is introduced to connect the potential function 
change across the interface. This method was first used by Vanden-Broeck [22], and later extended 
by Turner & Vanden-Broeck [6], Moni & King [11], Laget & Dias [23], Woolfenden & Parau [14] 
and others. 

The complex potential in each layer is represented by 



w ± {z) = (j) ± {z) + if{z) f 



(3.1) 



where 0 is the potential, ^ is the stream function and z = x + iy. We shall use a hodograph 
transformation exchanging the dependent and independent variables. Now, 0 ± are the 
independent variables along the interface and free surface. Without loss of generality, we choose 
x = Oat0 + = 0~=O and write 0 + =f(</>~). On the interface, we choose y\r = 0. It follows from 
the choice of the dimensionless variables that iff = —H on the rigid bottom and i/s = 1 on the free 
surface. We use 0 _ as the main variable. Therefore, the profiles of the interface and of the free 
surface are parametrized by £ _ (0 _ ), £ + (0 _ ), r]~((p~) and 1 + r) + {(j)~) r respectively. The lower layer 
is mapped into the infinite strip — H < ^ < 0 of the 0 _ + i^ plane. Using the method of images, 
the kinematic boundary condition on the horizontal bottom is satisfied by reflecting the lower 
layer into the bottom. This resulting lower layer is mapped into the infinite strip — 2H < ^ < 0 
of the 4>~ + ixj/ plane. We suppress the superscript in (p~ for the sake of simplicity. Applying the 
Cauchy integral formula to the lower fluid, we obtain 



L (<thr 



1 = 



2H&7 - 1) - (0 ■ 



■0^ 



(0"0m) 2 +4H 2 



d0 



-d0 



(3.2) 



For the upper layer, we need to apply the Cauchy integral formula to both the free surface and 
the interface. This yields the following two equations: 



tch It 



and 



1 = 



~U) + U(<t>) -f(4>m)]^ 



f(4>) -f(4>m) 



[/(</>) -f(4>m)] 2 + 1 



d<t> 



J(<t>)-f(cj> m ) 



(3-3) 



(3-4) 



The dynamic boundary conditions on the top free surface and the interface are used to close 
the system. In the transformed plane, the dynamic boundary condition takes the form 



where 

+ (it? 0-6) 



and 



Finally, the dynamic boundary condition on the interface takes the form 
F 2 



2 



o ^ o - (1 - R) 

&n 2 +te) 2 



+ (1 - R)^- = 0. (3.8) 



We seek solutions, symmetric about the centre of the wave, of the system (3.2)-(3.8) by a finite 
difference method. The numerical scheme for computing fully localized solitary waves is similar 
to [14], but modified for the present problem. The numerical computations are performed in a 
truncated long domain with length L. We introduce a uniform mesh 

0; = O'-l)^-^, i = l,2 N+l (3.9) 

and the corresponding unknowns 

% )<«- (!)<«■ 

Therefore, ^ and/^ are the functions of ry^, and, briefly speaking, they can be computed in the 
Newton iteration algorithm in the following order: 



computing %. via (3.2) 



computing via (3.8) 



computing via (3.4) 



We calculate ry^ , k and / by integrating their derivatives using the trapezoidal rule. The 
derivatives r\^, ^ and are computed by the five-point centred difference scheme. In fact, 
there are 3N unknowns, because (0n+i) = ffy (0n+i) = ^(0n+i) — 0 owing to the symmetry of 
the wave profiles. We then substitute all the expressions in (3.3), (3.5) and (3.7) evaluated at middle 
points 0 z m = {4>i + 0 z+ i)/2, i = 1,2, . . . ,N (a fourth-order interpolation formula is applied between 
the grid points and the middle points where necessary). This leads to 3N nonlinear algebraic 
equations for the 3N unknowns, which gives a closed system when the Froude number is 
given a priori. 

The nonlinear system is solved via a Newton iteration algorithm, and the program is 
considered to have converged to a solution when the L°° norm of the residual error is less than 
l(T 10 .The initial guess for Newton iteration is obtained by using the vanishing pressure method 
(see [24] for an example): an artificial fully localized pressure is applied either on the interface 
or on the free surface, then gradually eliminated using a continuation method by increasing the 
amplitude of the interface or free surface, and finally one obtains a solution without forcing. Once 
one solution is found, other solitary waves on the branch can be computed via a straightforward 
continuation method by changing the Froude number, density ratio, depth ratio or flexural 
rigidity as the parameter. Most of the computations were performed with the domain size L 
between 40 and 120, and the mesh size less than 0.1 to obtain sufficiently accurate numerical 
solutions. 

As Eb = 0, the governing equations (2.2)-(2.8) reduce to those of pure gravity interfacial waves 
with a free surface. In this case, barotropic elevation solitary waves have been numerically 
computed for the full Euler system by Moni & King [11] and Woolfenden & Parau [14]. To validate 
our numerical code, we checked the results produced by our program against the results of these 
authors. For R = 0.3, H = 1 and F = 1.404, Moni & King found the amplitudes of the interface and 
the free surface to be 0.134 and 1.552, whereas Woolfenden & Parau found them to be 0.137 and 
1.562, respectively. Our numerical results read 0.137 and 1.562, which are in excellent agreement 
with those of Woolfenden & Parau, confirming the accuracy of our numerical method. 



(a) (b) 




Figure 4. (0) Speed-amplitude bifurcation diagram for in-phase depression solitary waves with H = 1, E b = 1, R = 0.9 and 
the bifurcation point F = 1.2426 is predicted by linear theory (2.13); the upper branch is the free surface amplitude 1 + ^ + (0), 
and the lower branch is the interfacial amplitude rj~(0). (b) Typical wave profiles with F = 1.1473, 1 + ^ + (0) = 0.6 and 
r]~ (0) = —0.1776 (points labelled by squares in a); note that the free surface (a(\),b(\)) and the interface (000,600) are shown 
using different vertical scales. 

(ii) Fast (barotropic) solitary waves 

As shown in %2b, the dispersion curve of F + for the flexural-gravity problem always has a 
minimum at a non-zero wavenumber. Therefore, one can expect solitary waves with oscillatory 
decaying tails that bifurcate from infinitesimal periodic waves at the minimum of the external 
modes. It is noted that wavepacket barotropic solitary waves also appear in [14] when the term 
owing to elastic bending is replaced by surface tension. 

Two branches of solitary waves are found: depression waves with troughs at x = 0 and 
elevation waves with peaks at x = 0. Figures 4a and 5a show the bifurcation diagrams for these 
two branches with H = l, E\> = 1 and R = 0.9. The critical Froude number (i.e. the value of F 
corresponding to the minimum of the external mode) is then F + = 1.2426. Two typical profiles 
for F = 1.1473 and F = 1.138 are shown in figures 4b and 5b. They correspond to depression and 
elevation waves, respectively. At very small amplitudes, it is well known that waves bifurcating 
from extrema of a single mode (i.e. one free boundary) dispersion relation at finite k approach 
wavepackets of monochromatic waves of critical k with envelopes described by the modulational 
cubic nonlinear Schrodinger equation (NLS) [23]. Rigorous justifications of the NLS equations 
exist for single layer fluids [25]. When a free surface and an interface are both present, the 
NLS analysis is much more complicated and beyond the scope of this paper. By following the 
bifurcation branches to the highly nonlinear regime, we found that the waves ultimately approach 
configurations where the free surface and the interface tend to touch each other. Figure 6 shows 
typical profiles which have developed remarkable, almost touching structures for both depression 
(F = 0.2918) and elevation (F = 0.3) solitary waves. We also have conducted similar computations 
by fixing R = 0.9 and changing H and E\, by using a straightforward continuation method from the 
solutions shown. We found that the wave profiles and the bifurcation mechanism are qualitatively 
similar to those presented in figures 4-6. 

The limiting configurations of figure 6 differ from those found for one-layer gravity-capillary 
or gravity-flexural solitary waves, which ultimately converge to overhanging structures with 
intersection points for the depression branch [24,26] or exhibit snake-like bifurcation behaviour 
connecting to multi-packet wave profiles for the elevation branch [27,28]. 

(iii) Slow (baroclinic) solitary waves 

We now show that there is another type of solitary wave (baroclinic solitary waves) which 
bifurcate from the trivial uniform stream at values of the parameters corresponding to the maxima 




Figure 5. (a) Speed-amplitude bifurcation diagram for in-phase elevation solitary waves with H = 1, £ b = 1, R = 0.9 and 
the bifurcation point F = 1.2426 is predicted by linear theory (2.9); the upper branch is the free surface amplitude 1 + r] + (0) 
and the lower branch is the interfacial amplitude rj~(0). (b) Typical wave profiles with F = 1.138, 1 + r] + {0) = 1.1556 and 
rj~ (0) = 0.0518 (points labelled by squares in a); note that the free surface (tf(i),b(i)) and the interface (ff(ii),f>(ii)) are shown 
using different vertical scales. 



(*) (b) 




Figure 6. Highly nonlinear in-phase solitary wave profiles which are close to the extreme structure: (a) depression 
solitary waves (labelled by circles in figure 4a) with F = 0.2918, 1 + 77+ (0) = -0.81 and rj-(0) = -0.9002 and 
(b) elevation solitary waves (labelled by circles in figure 5a) with/ 7 = 0.3,1 + rj + (0) = 1.1927 and ^"(0) = -0.0308. 



in figure 2b. Unlike the external mode F + , there are two possibilities for the extreme value 
(global maximum) of F~: it can be attained at zero wavenumber or at non-zero wavenumber 
(as shown in figure 2b). The properties of the out-of-phase solitary waves for these two cases are 
significantly different. 

When the global maximum of F~ is attained at zero wavenumber, the bifurcation mechanism 
of the solitary waves is similar to those for interfacial long waves which, at small amplitude, 
are described by solutions of the Korteweg-de Vries equation. The waves are similar to the pure 
gravity interfacial solitary waves under the 'rigid lid' approximation [6,29], because the amplitude 
of the interface is relatively large compared with the displacement of the free surface. Figure 7 
presents bifurcation diagrams and typical profiles of the solitary waves with a set of parameters 
H=\,R = 0.9, E h = 0.5. These waves bifurcate from the uniform stream line at F = 0.1596 which 
coincides with the value predicted by the linear theory at k = 0, and the amplitude of the wave 
increases as the Froude number increases as shown in figure la. The waves are of baroclinic type 
(free surface and interface have different phases) with a positive amplitude on the interface. The 
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Figure 7. (a) Speed-amplitude bifurcation diagram for out-of-phase solitary waves with H = \, £b = 0.5, R = 0.9 and 
the bifurcation point F = 0.1596 is the minimum off - of the linear theory; the upper branch is the free surface amplitude 
1 + r] + {0) and the lower branch is the interfacial amplitude r]~{0). (b) The typical wave profiles and the broadening as 
the Froude number increases: F = 0.183594 (dotted line), F = 0.183606 (dashed line) and F = 0.183608 (solid line); the free 
surface profiles are shown at the top, whereas the interface profiles are shown at the bottom. 



most conspicuous feature is a broadening of the wave, namely the midsections of the interface 
and the free surface develop plateaus, when the Froude number approaches a limiting value. The 
plateaus can become infinitely long; therefore, the flow in the far field and the flow in the middle 
can be referred to as parallel conjugate flows [29,30]. From this point of view, we can calculate the 
limiting value of F by solving the following algebraic equation: 

(8 + R)F 6 - 12(1 + H)F 4 + 2(3 - R)(l + H) 2 F 2 - (1 - R)(l + H) 3 = 0. (3.11) 

This relation was first given by Dias & Il'ichev, and readers are referred to [16] for the detailed 
derivation (they derived (3.11) in the absence of an elastic cover, but their result remains valid 
for our problem, because the elastic sheet does not affect the result when the free surface and the 
interface are flat). The theoretic limiting value given by (3.11) is 0.1849, which is in agreement 
with the numerical result (~ 0.1836) as shown in figure 7a. It is worth mentioning that, in the 
broadening process, both the amplitude and the Froude number change slowly approaching the 
respective limiting values. This fact indicates that neither of them is an appropriate parameter in 
the continuation method, and, therefore, in numerical computations, we use the area beneath the 
interface as the variable parameter as suggested by Turner & Vanden-Broeck [6]. 

Long interfacial waves with a rigid lid may be of elevation or depression depending on 
whether the equilibrium interface is above or below a particular depth (the midline in the limit 
as the densities are close; [29]). Thus, when we increase the relative depth H, the interfacial 
amplitude becomes negative, whereas the free surface and the interface are still in opposite 
phases. In figure 8a, we pick H = 3, R = 0.9 and E\> = 0.1, and plot the speed-amplitude bifurcation 
diagrams. The out-of-phase solitary waves under this set of parameters also display the 
broadening phenomenon (see the typical wave profiles in figure 8b). The numerical limiting value 
of F is about 0.3222, which coincides with the theoretical prediction (^ 0.3203) calculated by (3.11). 

We note that the elastic bending has little effect in the computation of the long solitary waves 
given the small variation of the free surface: we tested the profiles of the solitary waves under the 
set of parameters H = 3, R = 0.75 and F = 0.51 by varying E\>. It turns out that the wave profiles 
do not change appreciably, with the absolute difference between two profiles computed with 
different being less than 10 -3 for E^ varying in the range [0.01,1]. This is to be contrasted 
with the long solitary waves of the external mode [12,13]. There, because the elastic bending 
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Figure 8. (a) Speed-amplitude bifurcation diagram for out-of-phase solitary waves with H = 3, £b = 0.1, R = 0.9 and 
the bifurcation point F = 0.277 is the minimum off - calculated by (2.14); the upper branch is the free surface amplitude 
1 + ?7 + (0), and the lower branch is the interfacial amplitude rj~(0). (b) The broadening of the wave profile as the Froude 
number increases: F = 0.3215 (dotted line),/ 7 = 0.322140 (dashed line) and F = 0.322204 (solid line); the free surface (a(\),b(\)) 
and the interface (<7(ii),6(ii)) are shown separately using different vertical scales. 



dramatically changes the dispersion curve of the external mode for large wavenumbers, non- 
decaying periodic tails are ruled out in the far field of long waves. 

For the case when the maximum of the internal mode appears at a non-zero wavenumber, 
the bifurcation mechanism and the wave profiles are completely different from those mentioned 
above, in this section, and are similar to those mentioned in the previous section bifurcating 
from the minimum of the external mode. Consider the numerical experiment in figure 9. This 
experiment is carried out with H = 3, R = 0.1 and E\, = 2, and therefore the maximum 0.9451 of F~ 
occurs atk = 0.7052. Typical wave profiles are shown in figure 9b which are wavepacket solitary 
waves featuring damped oscillatory tails. From the speed-amplitude bifurcation diagrams 
presented in figure 9a, we observe that the amplitude of this type of solitary wave increases 
as the Froude number increases, and that the branch fills the whole gap of speeds between the 
maximum of the internal mode and the minimum of the external mode. These solitary waves 
do not broaden, and thus the theoretical prediction of the limiting value of the Froude number 
(1.1547 under those parameters) is not applicable. It is to be noted that no in-phase solitary waves, 
which would bifurcate from infinitesimal periodic waves at the minimum of F + , are found in our 
numerical experiments under these parameters. 

(b) Generalized solitary waves 
(i) Series truncation method 

Numerically generalized solitary waves are computed via the modified version of the method 
stated in [13]. We approximate a generalized solitary wave by long periodic waves with period 
L. As done in §3a(i), the fluid geometry is mapped to the 0 ± — y\r domain with the upper layer 
[-L/2,L/2] x [0,Q + ] and the lower layer [-L/2,L/2] x [0,-Q - ], where Q ± are two unknown 
constants corresponding to the values of x/r on the free surface and the flat bottom. All the 
solutions can be approximated using truncated Fourier series, namely in the lower layer 
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Figure 9. (a) Speed-amplitude bifurcation diagram for out-of-phase solitary waves with H = 3, ft = 2, R = 0.1 together 
with two bifurcation points F = 0.9451 — the maximum off - — and F = 1.2980 — the minimum of F + ; the upper branch is 
the free surface amplitude 1 + r] + {0), and the lower branch is the interfacial amplitude ^~(0). (b) The typical wave profile 
which is labelled by circles in (a): F = 1.2477, 1 + r? + (0) = 1.56 and rj~{0) = -0.3257. 
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Because there is a potential jump across the interface, we express 0 ± as 
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Therefore, all the unknowns are parametrized by the variable £. It is obvious that, as iff = — Q , 
one has y~ = —Q~ + ao = —H, which gives Q~ . Overall, we have the following 4N + 4 unknowns: 
a n (N + 1), b n (N + 1), c n (N - 1), d n (N), B ± and Q+. Accordingly, using the mesh (3.9), one can have 
N — 1 equations from x~ = x + at the interface, N + 1 equations from y~ = y + at the interface, 
2N + 2 equations from the dynamic boundary conditions on the interface and free surface. 
Additionally, two mean-value conditions complete the problem, 
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We should remark that this series truncation method is not adequate to compute large-amplitude 
waves for our problem. As pointed out in [13], the computation of large-amplitude waves requires 
a large amount of mesh points, and, using this numerical scheme, the terms cosh(27rn^/L) and 
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Figure 10. (a) Typical profiles of long periodic waves with non-decaying periodic tails with parameters R = 0.1, H = 3,E b = 
0.02 and F = 1; superposition of the waves computed in different domains (L/2 = 70 for the solid line and L/2 = 76.3 for 
the dashed line) suggests the existence of the generalized solitary waves as L -> oo. (b) Profiles of long free-surface and 
interfacial internal waves for R = 0.1, E b = 0.1, H = 3, L/2 = 70 and three values of the Froude number: F = 0.95 (dotted 
line), 1.00 (solid line), 1.03 (dashed line), (c) Profiles of generalized solitary waves for R = 0.1, E b = 0.1, F = 1, L/2 = 70 and 
two different depth ratios: H = 0.22 (solid line), H = 0.28 (dashed line), (d) Profiles of the free surface and the interface for 
R = 0.1, H = 3, F = 1, L/2 = 70 and two flexural— rigidity constants: £ b = 0.11 (solid line), £ b = 0.033 (dashed line). 



sinh(27rn^/L) in x + and y + grow exponentially as n increases. Thus, in §3fr(ii), we focus mainly 
on relatively small-amplitude generalized solitary waves, and their complicated tails owing to 
the effect of elastic bending. 

(ii) Some numerical results 

A generalized solitary wave consists of a localized midsection and non-decaying oscillatory tails. 
This type of wave arises from the scenario that a long wave [12,13,17,31] or a finite-amplitude 
wavepacket solitary wave [32] resonates with other linear waves which propagate at the same 
speed. Here, we consider the parameter regime in which long waves bifurcating from F~ at k = 0 
resonate with finite wavenumber waves from the external branch near F + . In order for this to 
happen before the internal waves broaden and reach a limiting amplitude, the gap between the 
maximum of F~ and the minimum of F + must be sufficiently small (or non-existent). We note that 
this requires parameters that are unphysical for oceanic applications. 



Generalized solitary waves, which mathematically extend to infinity, are approximated by long 
periodic waves in the numerical experiments. As an example, we fix R = 0.1 and first compute 
approximate long solitary waves as periodic waves with long flat troughs. Then, we vary the 
parameters H, F and to modify the dispersion relation until a resonance occurs resulting in the 
appearance of periodic tails. For example, we first compute the wave with R = 0.1, H = 3, E^ = 0.1, 
L/2 = 70 and choose the Froude number F = 0.95 which is in the gap between the minimum of 
F + and the maximum of F~ in the dispersion relation. The dotted line in figure 10b shows that, 
for this Froude number, the wave is solitary-like and flat in the far field. We then increase the 
Froude number and ripples appear on the tails of the profiles as expected (see the dashed and 
solid lines in figure 10b) when F exceeds the minimum of F + . 

In figure 10a, we present typical free-surface profiles and corresponding interfaces for the 
parameters R = 0.1, H = 3, F = 1, E\> = 0.02, and the superposition of waves computed within two 
different domains. The midsection of the free surface and that of the interface are in opposite 
phases. It shows that, as the domain size increases, more and more periodic waves can be added 
in the far field without a notable change of the portion of the profile already obtained for smaller 
domain sizes. The computation strongly suggests the existence of generalized solitary waves in 
the limit as the size of the domain approaches infinity. Figure lOd demonstrates the formation 
of a trough at the origin of the free surface as the flexural-rigidity constant varies, indicating 
that curvature at the origin may change sign along with the changing parameters. There are 
many different families of generalized solitary waves characterized by an infinitely long train of 
ripples in the far field. The periodic tails of generalized solitary waves in a two-layer system with 
elastic cover are more complicated than that of the waves without an elastic sheet (i.e. E^ = 0). 
The solid lines in figure 10c ,d show a Wilton-ripple-like periodic tail rather than the sinusoid- 
like tail of the pure gravity case. That is because F + is not monotonic, and the finite-amplitude 
solitary waves can resonate with two linear periodic waves when the Froude number is above 
the minimum of F~ . Generalized solitary waves with Wilton ripple far fields do not appear to 
have been calculated before. We conjecture that a similar phenomenon can also be observed in 
a two-layer system with a free surface, when the surface tension or interfacial tension is taken 
into account. 



4. Conclusions 

In a two-layer fluid system with a free upper surface, there are two distinct dispersion curves, 
F + for the external mode and F~ for the internal mode. Depending on the specific physical 
effects, either on the free surface or on the interface, there may be a gap between F + and 
F~ where varieties of gap solitary waves can exist. In this paper, we considered interfacial 
gravity waves with an elastic cover on the top of the upper layer fluid — a problem of interest 
in oceanic applications [4]. The nonlinear elastic model is based on the special Cosserat theory 
of hyperelastic shells. Steady, symmetric gap solitary waves are computed via a boundary 
integral method. For the realistic density ratios, there are two types of fully localized travelling 
waves: barotropic and baroclinic solitary waves. Barotropic solitary waves behave like one-layer 
free surface waves with wavepacket solitary waves, including depression and elevation waves, 
existing in this case. Both these are found to approach a limiting configuration that the free surface 
and the interface tend to touch each other. Baroclinic solitary waves, bifurcating from infinitesimal 
internal long periodic waves at Froude number F~ (0), are restricted by a critical Froude number 
corresponding to parallel conjugate flows. If this critical Froude number is between the minimum 
of F + and F~(0), the midsections of the interface and the free surface develop plateaus which 
become wider as the Froude number approaches its critical value from below. Baroclinic wave- 
packet solitary waves, bifurcating from finite wavenumber, are found for smaller density ratios 
when the maximum of F~ may occur at finite wavenumber. These solitary waves with decaying 
oscillatory tails are found numerically to be no longer confined by the critical Froude number, 
and can fill the entire gap between the minimum of F + and the maximum of F~ . 



Small-amplitude generalized solitary waves are computed using a series truncation method 
and for the small density ratio. The Wilton-ripple-like periodic tails owing to the non- 
monotonicity of the external mode F + , which are distinct from the case when flexural effects 
are ignored, are highlighted. While it is beyond the scope of this paper, the computation of large- 
amplitude generalized solitary waves for the irrotational Euler system demands more refined 
numerical schemes and could be particularly interesting. 

The stability of steady interfacial waves with free surface is an open problem and merits further 
investigation. The problem is complicated by the fact that interfacial waves induce a jump in 
tangential velocity across the interface resulting in small-scale shear instabilities. Despite this, 
we expect that waves similar to those we computed could be found in oceanic conditions or 
laboratory experiments as is the case of interfacial waves with a free surface in the presence of 
gravity alone. Furthermore, the dead-water phenomenon is likely to be observed when the load 
on the elastic sheet moves slowly enough to excite the internal modes. 
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